A closer look at time averages of the logistic map at the edge of chaos 
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The probability distribution of sums of iterates of the logistic map at the edge of chaos has been 
recently shown [see U. Tirnakli, C. Beck and C. Tsallis, Phys. Rev. E 75, 040106(R) (2007)] to 
be numerically consistent with a g-Gaussian, the distribution which, under appropriate constraints, 
maximizes the nonadditive entropy S q , the basis of nonextensive statistical mechanics. This analysis 
was based on a study of the tails of the distribution. We now check the entire distribution, in 
particular its central part. This is important in view of a recent (/-generalization of the Central 
Limit Theorem, which states that for certain classes of strongly correlated random variables the 
rescaled sum approaches a g-Gaussian limit distribution. We numerically investigate for the logistic 
map with a parameter in a small vicinity of the critical point under which conditions there is 
convergence to a q-Gaussian both in the central region and in the tail region, and find a scaling 
law involving the Feigenbaum constant 8. Our results are consistent with a large number of already 
available analytical and numerical evidences that the edge of chaos is well described in terms of the 
entropy S q and its associated concepts. 

PACS numbers: 05.20.-y, 05. 45. Ac, 05.45.Pq 

One of the cornerstones of statistical mechanics and of probability theory is the Central Limit Theorem (CLT). 
It states that the sum of N independent identically distributed random variables, after appropriate centering and 
rescaling, approaches a Gaussian distribution as N —> oo. In general, this concept lies at the very heart of the 
fact that many stochastic processes in nature which consist of a sum of many independent or nearly independent 
variables converge to a Gaussian process [l], 0] • On the other hand, there are also many other occasions in nature for 
which the limit distribution is not a Gaussian. The common ingredient for such systems is the existence of strong 
correlations between the random variables, which prevent the limit distribution of the system to end up being a 
Gaussian. Recently, for certain classes of strong correlations of this kind, it has been proved that the distribution of 
the rescaled sum approaches a q-Gaussian, which constitutes a ^-generalization of the standard CLT This 
represents a progress since the q-Gaussians are the distributions that optimize the nonadditive entropy S q (defined to 
be S q = (1 — J2iPl) I il ~~ 1))' 011 which nonextensive statistical mechanics is based A ^-generalized CLT was 

expected for several years since the role of q-Gaussians in nonextensive statistical mechanics is pretty much the same 
as that of Gaussians in Boltzmann-Gibbs statistical mechanics. Therefore it is not surprising at all to see q-Gaussians 
replace the usual Gaussian distributions for those systems whose agents exhibit certain types of strong correlations. 

Immediately after these achievements, an increasing interest developed for checking these ideas and findings in real 
and model systems whose dynamical properties make them appropriate candidates to be analyzed along these lines. 
Cortines and Riera have analysed stock market index changes for a considerable range of time delays using Brazilian 
financial data @ and found that the histograms can be well ap prox imated by q-Gaussians with q m 1.75 (see Fig. 6 of 
0). Another interesting study has been done by Caruso et al. l(j] by using real earthquake data from the World and 
Northern California catalogs, where they observed that the probability density of energy differences of subsequent 
earthquakes can also be well fitted by a q-Gaussian with roughly the same value of q, i.e., q w 1.75 (see Fig. 2 of flOj] ^> . 
A more recent contribution along these lines consists in a molecular dynamical test of the q-CLT in a long-range- 
interacting many-body classical Hamiltonian system known as HMF model [ill ], where it was numerically shown that, 
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in the longstanding quasi-stationary regime (where the system is only weakly chaotic) , the relevant densities appear 
to converge to g-Gaussians with q ss 1.5 (see for example Fig. 7 of fill, see also [Hj]). Moreover, q-Gaussians have 
also been observed in the motion of Hydra cells in cellular aggregates |13j ] , for defect turbulence [14| , silo drainage of 
granular matter [l5j], cold atoms in dissipative optical lattices [l6|, and dissipative 2D dusty plasma [T^. Finally, in 
a recent paper, we numerically investigated the central limit behavior of deterministic dynamical systems (l8j . where 
one of our main purposes was to see what kind of limit distributions emerge for the attractor whenever the dynamical 
system is not mixing (for example at the edge of chaos, where the Lyapunov exponent vanishes) and thus the standard 
CLT is not valid anymore. In 18], using the well-known standard example of discrete one-dimensional dissipative 
dynamical systems, the logistic map, defined as 

x i+ i = l-ax\ , (1) 

(where < a < 2; |x t | < 1; t = 0, 1, 2, ... ), we numerically checked that, at the edge of chaos (i.e., close to the critical 
parameter value a c = 1.401155189092...), the tails of the limit distribution were consistent with a g-Gaussian having, 
once again, a value of q close to 1.75. However, the central part of the distribution was not meticulously studied, 
and neither was studied the precise dependence on the distance a — a c and on the iteration number. In the present 
manuscript, our aim is to focus on these points, having a closer look at sums of iterates of the logistic map close to 
its chaos threshold. 

Although the iterates of a deterministic dynamical system can never be completely independent, one can still 
prove some standard CLTs for such systems [ill US IH|: provided that the assumption of independent identically 
distributed random variables is replaced by the property that the system is sufficiently mixing (i.e., asymptotic 
statistical independence). As an example one can consider the logistic map at a = 2 where it is strongly mixing. For 
this system, it can be rigorously proved [19|, l2fj that the distribution of the quantity 

JV 

y:=5>i-<z» (2) 

i=l 

becomes Gaussian for N — > oo after appropriate rescaling with a factor 1/yN, regarding the initial value x\ as a 
random variable with a smooth probability distribution. Here (x) denotes the mean of x, which happens to vanish 
for the special case a — 2. This is a highly nontrivial result since the iterates of the logistic map at a — 2 are not 
independent but exhibit complicated higher-order correlations described by forests of binary trees [22J . Gaussian limit 
behavior is also numerically observed for other typical parameter values in the chaotic region of the logistic map [1 81 ] . 
Indeed, whenever the Lyapunov exponent of the one-dimensional map is positive, one expects the CLT to be valid 

a 

Now we are ready to discuss the behavior of the logistic map at the edge of chaos for which a standard CLT is not 
valid due to the lack of mixing. In order to calculate the average in Eq. ((2]), it is necessary to take the average over a 
large number of N iterates as well as a large number of randomly chosen initial values , namely, 

2=1 i=l 

These conditions are important due to potentially non-ergodic and non-mixing behavior. 

In principle, at the edge of chaos, taking N — > oo is not the only ingredient for the system to attain its limit 
distribution. It is also necessary to localize the critical point (the chaos threshold) with infinite precision. In other 
words, theoretically, for a full description of the shape of the distribution function on the attractor, one needs to take 
the a c value with infinite precision as well as taking N — > oo. On the other hand, in numerical experiments, neither 
the precision of a c nor the N values can approach infinity. In fact, numerically one can see the situation as a kind of 
interplay between the precision of a c and the number N of iterates. For a given finite precision of a c (slightly above 
the exact critical value), if we use a very large N, then the system quickly feels that it is not exactly at the chaos 
threshold, and the central part of its probability distribution function typically becomes a Gaussian (with only small 
deviations in the tails). On the other hand, for the same distance to a c , if we take N too small, then the summation 
given by Eq. ([2]) starts to be inadequate to approach the edge-of-chaos limiting distribution, and the central part of 
the distribution around zero exhibits peaks. This is indeed a direct consequence of the fact that the attractor of the 
system at the edge of chaos is a fractal that only occupies a tiny part of the full phase space (see [23| for details). 
This is the reason why the central parts of the distributions shown in [l8| do not present the typical smooth shape of 
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FIG. 1: Probability density function P(y) of the quantity y rescaled by -P(O). The map is close to the edge of chaos with 5 
digits precision (a = 1.40116). The tendency to approach the Gaussian is evident as N increases and a is kept fixed. 

g-Gaussians. Indeed, the values of N chosen in [lj| (2 14 and 2 15 ) are too small for the precision of a c (1.401155189092) 
to obtain a complete picture of the entire distribution including both central parts and tails. On the other hand, 
for the above precision of a c , one can think about numerical experiments with N values at the level of, say, 2 or 
more, for which the central part would approach a Gaussian since the system starts to realize that it is not exactly 
at the edge of chaos. We observe that between these two extremes, for a given precision of a c , there exists a range of 
values of N for which the probability density of the system is well approximated by a ^-Gaussian in the entire region. 
Unfortunately, these kinds of large N values which are necessary to fully verify this observation if we approach a c 
with say 12 digits precision (as we did in [l8l |) cannot normally be reached in numerical experiments. However, we can 
check this scenario using less precision for a c . This will in turn make the appropriate N value become small enough 
so that we are able to handle the numerics with standard computers. 

As a representative illustration, we first focus on an a value in the vicinity of the critical point with 5 digits precision 
(a = 1.40116). For this case (and several other cases which are not described here), we numerically verified the above- 
mentioned scenario, as can be seen in Fig. 1. In our simulations, after omitting a transient of the first 2 12 iterates (we 
checked that the results are independent of the omitted transient length as long as it is large enough), we calculated 
the quantity y in Eq. (|2|) for various N values and obtained its probability distribution from an ensemble of uniformly 
distributed initial values. It is worth mentioning that a similar picture emerges for almost any vicinity of the critical 
point, but of course with different N values. It is clearly seen from the figure that, for a = 1.40116, N — 2 19 is so 
large that the density approaches a Gaussian in the central region with small deviations in the tails (further increase 
of N values would make the whole curve become a Gaussian), whereas for TV = 2 12 the curve has heavy tails and a 
peaked central part (this curve can be fitted neither by a Gaussian nor by a g-Gaussian in the entire region). On the 
other hand, between these two extreme cases, there is an appropriate range of N values (around 2 17 for this example 
of a — a c ) for which the distribution is consistent with a g-Gaussian of the form 
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(where q and (3 are suitable parameters) for the entire region. 

Let us now provide a theoretical argument what the optimum value of N could be to achieve best convergence to 
a g-Gaussian. Finite precision of a c means that the parameter a of the system is at some distance \a — a c \ from the 
exact critical point a c — 1.401155189092.... Suppose we are slightly above the critical point (a > a c ), by an amount 



l a ~ ac I^^T' ( 5 ) 

where 5 = 4.6692011... is the Feigenbaum constant. Then there exist 2™ chaotic bands of the attractor with a 
selfsimilar structure, which approach the Feigenbaum attractor for n — > oo by the band splitting procedure (see e.g. 
[24| . p. 10, for more details). Suppose we perform 2™ iterations of the map for a given initial value with a parameter 
a as given by Eq. ([5]). Then after 2™ iterations we are basically back to the starting value, because we fall into the 
same band of the band splitting structure. This means the sum of the iterates X)i=i x i essentially approach a 
fixed value w = 2 n (x) plus a small correction Aiui which describes the small fluctuations of the position of the 2 n th 
iterate within the chaotic band. Hence 

2" 

yi=^2(xi- (x)) = Awi, (6) 
If we continue to iterate for another 2 n times, we obtain 

2" +1 

y 2 = ^2 (ii - (x)) = Aw 2 . (7) 

i=2"+l 

The new fluctuation Aw2 is not expected to be independent from the old one Aw± , since correlations of iterates decay 
very slowly if we are close to the critical point. Continuing, we finally obtain 

4 n 

i=4"-2" + l 

if we iterate the map 4™ times in total. The total sum of iterates 



4" 2" 
i=l 3=1 

can thus be regarded as a sum of 2 n strongly correlated random variables Auij , each being influenced by the structure 
of the 2" chaotic bands at distance a— a c ~ 5~ n from the Feigenbaum attractor. There is a 1-1 correspondence between 
these 2™ random variables Awj,j = 1, . . . , 2" and the 2" chaotic bands of the attractor, which remains preserved if 
n is further increased. It is now most reasonable to assume that the above system of a sum of 2™ correlated random 
variables Awj exhibits data collapse (and hence convergence to a well-defined limit distribution) under successive 
renormalization transformations n^n+l^n + 2^n + 3- -- . The limit distribution may indeed be a q-Gaussian, 
as indicated by our numerical experiments. The above scaling argument implies that the optimum iteration time N* 
to observe convergence to a g-Gaussian limit distribution is given by 



N* - 2 2n (10) 
where, at a given distance a — a c , the number n is given by 



log I a — a c I . , 

n« x ■ n 

log S 

Another way to formulate our scaling argument is as follows. Consider a given distance from the critical point 
where there are 2 n chaotic bands. The relevant function in the Feigenbaum renormalization scheme is the 2™-th 
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TABLE I: Various a values used in the simulations and the associated values of n and TV*. 



a 


\d — &c | 


n 


N* 


1.40209 


9.348 • 10~ 4 


4.526~ 9/2 


2 y 


1.401588 


4.328 ■ 10~ 4 


5.026~ 10/2 


2 iu 


1.401354 


1.988 ■ 10" 4 


5.531~ 11/2 


2 11 


1.401248 


9.281 • 10~ b 


6.025~ 12/2 




1.401198 


4.281 • 10 _a 


6.527~ 13/2 


2 u 


1.401175 


1.981 ■ 10 _b 


7.027~ 14/2 


2 14 


1.4011644 


9.211 ■ 10~ B 


7.524~ 15/2 


2 lb 


1.40115945 


4.261 ■ HT° 


8.025~ 16/2 


2 lb 


1.40115716 


1.971 ■ 10 _b 


8.525~ 17/2 
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iterated function /* := f 2 rather than the original function / itself. The iterates of /* will be highly correlated. 
Let k be the number of iterates of /* that are being added up. For k >> 2 n we get a Gaussian distribution for the 
probability distribution of the sum, since the system feels that it is in the chaotic regime a > a c . For k << 2™ there is 
no chance to get anything smooth for the probability distribution of the sum since the iteration number is too small. 
The interesting intermediate case is the case k = 2™. Here, due to the strong correlations, there is the chance to get 
a smooth g-Gaussian limit distribution. Since in this case the number of chaotic bands is the same as the number 
of iterates of that are being added up, the theory is invariant under successive renormalization transformations 
n — *■ n. + 1 — »n. + 2.... But iteration number k — 2™ for /» corresponds to iteration number N* — 2 2n for the original 
map, which is our equation (JTUJ). 

In order to check these types of scaling arguments, we numerically studied various a values as listed in the Table. 
First we show in Fig. 2a and 2b the densities obtained for the case N = N* = 2 2n for each value of a, where 2n is 
odd and even, respectively. Three important aspects are evident: 

(i) the curves obtained for these cases exhibit a very clear data collapse; 

(ii) the envelopes of the histogram data can be well fitted everywhere (i.e. both in the tails and in the central region) 
by a g-Gaussian given by Eq.(0| with q — 1.68 if 2n is odd and q — 1.70 if 2n is even, for the respective a values 
given in the Table. It is also evident from this figure that small log-periodic modulations are present on top of the 
curves. The existence of such modulations is expected due to discrete scale invariance in these types of systems, as 
demonstrated long ago by de Moura et al. [25| (see also [26j]); 

(iii) as we approach a c more closely, and consequently the appropriate value N* increases, the region consistent with 
the g-Gaussian grows in size. In order to illustrate this statement we include Fig. 3 which shows the same data in a 
different representation. One can better see there how the g-Gaussians with log-periodic-like modulations develop as 
a approaches a c . g-Gaussians correspond to a straight line in these types of plots. We also determined and show in 
Fig. 4 how P(0) evolves with N for the cases studied in Fig. 2. 

In Fig. 5 we plot the g-logarithm (defined to be the inverse function of the g-exponential given in Eq. namely 
ln q (x) — (x 1 ^ 9 — 1)/(1 — q)) of the same data as in Fig. 2a. This provides a further visualization of the aspect 
already mentioned in (iii), namely that as the precision of a c and the value of N increase, the region consistent with 
a g-Gaussian extends in size. There is a clear numerical indication that g-Gaussians are a good approximation of the 
data if both the precision of a c and the value of N go to infinity. 

We also checked that our results are not induced by roundoff or similar numerical artifacts. Throughout our 
simulations we used double precision of Intel Fortran, achieving good statistics {rii n i — 10 8 for all cases). We also 
tested our results with less statistics but using higher precision (quadruple precision of Intel Fortran). Since no 
significant differences were observed, we used double precision with high statistics in most of our simulations, which 
allowed us to see the observed g-Gaussians to be modulated by small log-periodic-like oscillations. The detailed 
structure of these oscillations depends on the precise value of a in a very complicated way, which is to be expected, 
since the logistic map is well-known to exhibit very complex behavior as a function of a. Nevertheless, our numerical 
experiments provide evidence that the envelope of the data is always very well approximated by a q-Gaussian, provided 
the typical number of iterations is chosen according to the scaling relations given by Eqs. (fT"0|) and (fTTj) . Another 
interesting subject is, no doubt, to analyse the cases N « N* and N » N* in order to understand the crossover 
phenomena from the peaked region to g-Gaussian region and finally to the normal Gaussian region. This will be 
addressed elsewhere in the near future. 

Summarizing, we have presented numerical evidence that the distributions of sums of iterates of the logistic map 
in a close vicinity of the edge of chaos are well approximated by g-Gaussian probability distributions, provided the 
typical number of iterations scale in line with Eqs. (fTU|) and (fTTj) . This illustrates the strongly correlated nature of 
this paradigmatic nonlinear dynamical system (which models a great variety of more complex physical situations, 
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FIG. 2: Data collapse of probability density functions for the cases N = 2 2n , where 2n is odd (a) and even (b). As n increases, 
a good fit using a g-Gaussian with q — 1.68 and j3 = 6.2 (a) and q = 1.70 and f3 — 6.2 (b) is obtained for regions of increasing 
size. Inset: The linear-linear plot of the data for a better visualization of the central part. 



as it is well-known in the literature). The q-Gaussians are precisely the limit distributions of an important class 
of strongly correlated random systems in the realm of the recently proved g-generalized Central Limit Theorem. 
This feature, together with the fact that they optimize within appropriate constraints the nonadditive entropy S q , 
is of interest for the mathematical foundations of nonextensive statistical mechanics, as well as for many real-world 
problems that consist of sums of correlated random variables. Needless to say that the analytical proof of the results 
numerically obtained here would be most welcome. As open questions we may mention (i) the careful study of other 
dissipative and conservative maps, starting with the one-dimensional z-logistic one, and (ii) the investigation of the 
precise convergence radius to the g-Gaussian limit form and its oscillating small correction terms as a function of N 
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FIG. 3: Probability density functions plotted against 1 + (q — l)(3[yP(0)] 2 on a log-log plot for the cases N — 2 2n , where 2n is 
odd (a) and even (b). A straight line is expected with a slope 1/(1 — q) if the curve is a g-Gaussian. It is clearly seen how the 
straight line is surrounded by the log-periodically modulated curves. 



and a — a c . A full numerical study of these points certainly requires high computational power. 
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FIG. 4: iV-dependence of the rescaling factor P(0) that yields data collapse as displayed in Fig. 2a and 2b. 
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